Haemosphere

library(edgeR)
library(limma)

library(SingleCellExperiment)
library(SingleR)
library(Seurat)


library(biomaRt)

library(dplyr)
library(reshape2)

library(patchwork)
library(ComplexHeatmap)
library(ggplotify)
library(gridExtra)
library(grid)
library(viridis)

Loading required package: limma

Loading required package: SummarizedExperiment

Loading required package: MatrixGenerics

Loading required package: matrixStats


Attaching package: ‘MatrixGenerics’


The following objects are masked from ‘package:matrixStats’:

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
    colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
    colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
    colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
    colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
    colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
    colWeightedMeans, colWeightedMedians, colWeightedSds,
    colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
    rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
    rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
    rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
    rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
    rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
    rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
    rowWeightedSds, rowWeightedVars


Loading required package: GenomicRanges

Loading required package: stats4

Loading required package: BiocGenerics

Loading required package: parallel


Attaching package: ‘BiocGenerics’


The following objects are masked from ‘package:parallel’:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB


The following object is masked from ‘package:limma’:

    plotMA


The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs


The following objects are masked from ‘package:base’:

    anyDuplicated, append, as.data.frame, basename, cbind, colnames,
    dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
    grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
    rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which.max, which.min


Loading required package: S4Vectors


Attaching package: ‘S4Vectors’


The following objects are masked from ‘package:base’:

    expand.grid, I, unname


Loading required package: IRanges

Loading required package: GenomeInfoDb

Loading required package: Biobase

Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.



Attaching package: ‘Biobase’


The following object is masked from ‘package:MatrixGenerics’:

    rowMedians


The following objects are masked from ‘package:matrixStats’:

    anyMissing, rowMedians



Attaching package: ‘SingleCellExperiment’


The following object is masked from ‘package:edgeR’:

    cpm


Attaching SeuratObject


Attaching package: ‘Seurat’


The following object is masked from ‘package:SummarizedExperiment’:

    Assays



Attaching package: ‘dplyr’


The following object is masked from ‘package:biomaRt’:

    select


The following object is masked from ‘package:Biobase’:

    combine


The following objects are masked from ‘package:GenomicRanges’:

    intersect, setdiff, union


The following object is masked from ‘package:GenomeInfoDb’:

    intersect


The following objects are masked from ‘package:IRanges’:

    collapse, desc, intersect, setdiff, slice, union


The following objects are masked from ‘package:S4Vectors’:

    first, intersect, rename, setdiff, setequal, union


The following objects are masked from ‘package:BiocGenerics’:

    combine, intersect, setdiff, union


The following object is masked from ‘package:matrixStats’:

    count


The following objects are masked from ‘package:stats’:

    filter, lag


The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union


Loading required package: grid

========================================
ComplexHeatmap version 2.8.0
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference

If you use it in published research, please cite:
Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
  genomic data. Bioinformatics 2016.

The new InteractiveComplexHeatmap package can directly export static 
complex heatmaps into an interactive Shiny app with zero effort. Have a try!

This message can be suppressed by:
  suppressPackageStartupMessages(library(ComplexHeatmap))
========================================



Attaching package: ‘gridExtra’


The following object is masked from ‘package:dplyr’:

    combine


The following object is masked from ‘package:Biobase’:

    combine


The following object is masked from ‘package:BiocGenerics’:

    combine


Loading required package: viridisLite

# Set working directory to project root
setwd("/research/peer/fdeckert/FD20200109SPLENO")

# Source files
source("plotting_global.R")

Loading required package: RColorBrewer

Parameter

# Seurat files 
cnt_file <- "data/haemosphere/Haemopedia-Mouse-RNASeq_raw.txt"
meta_file <- "data/haemosphere/Haemopedia-Mouse-RNASeq_samples.txt"

so_file <- "data/object/seurat.rds"

filter_min_cpm = 0.5
filter_min_expressed_samples = 2 
normalization_method = "TMM"
adj_p_cutoff = 0.05
min_rows = 50

# Plotting Theme
ggplot2::theme_set(theme_global_set()) # From project global source()

# Settings 
options(digits = 4)

Loading required package: ggplot2

Helper functions

# Takes a prcomp pca object and group labels as input to compute a PCA scattern plot of the first two PCA 
pca_scattern <- function(pca, group, group_color) {
    
    pca_data <- cbind(pca$x, data.frame(group = group))
    pca_scattern <- ggplot(pca_data, aes(x = PC1, y = PC2, fill = group)) + 
      geom_point(color = "black", pch = 21, size = 3) + 
      ggtitle("log2(CPM) - z-score") + 
      xlab(paste("PC1 (", round(100 * (pca$sdev^2)[1] / sum(pca$sdev^2), digits = 1), " %)")) + 
      ylab(paste("PC2 (", round(100 * (pca$sdev^2)[2] / sum(pca$sdev^2), digits = 1), " %)")) +
      scale_fill_manual(values = group_color) +
      guides(fill = guide_legend(title = "Group")) + 
      theme(aspect.ratio = 1, 
            panel.border = element_blank(), 
            axis.line    = element_line())
    
    return(pca_scattern)
}

top_hm <- function(deg, cnt = NULL, meta = NULL, level = "", contrast = "", top = 50, min_contrast_log2cpm = 2) {

    # Filter and sort cnt by meta samples 
    # Columns of cnt must be corespondend to meta sample id to compute mean cpm by group
    cnt <- cnt[, colnames(cnt) %in% meta$sample_id]  
    cnt <- cnt[, match(meta$sample_id, colnames(cnt))]
    
    # Compute cpm and log2(cpm)
    cpm <- edgeR::cpm(cnt, normalized.lib.sizes = FALSE, log = FALSE, prior.count = 0)
    log_cpm <- edgeR::cpm(cnt, normalized.lib.sizes = FALSE, log = TRUE, prior.count = 1)
    
    # Define grouping levels to compute the group mean log2 cpm
    if(level == "main_labels") {group <- meta$main_labels}
    if(level == "fine_labels") {group <- meta$fine_labels}
    
    # Mean group cpm and log2(cpm)
    colnames(cpm) <- group    
    cpm <- lapply(unique(group), function(i) data.frame(X = rowMeans(cpm[, grepl(i, colnames(cpm), fixed = TRUE)])))
    cpm <- do.call("cbind", cpm)
    colnames(cpm) <- unique(group)
                  
    colnames(log_cpm) <- group
    log_cpm <- lapply(unique(group), function(i) data.frame(X = rowMeans(log_cpm[, grepl(i, colnames(log_cpm), fixed = TRUE)])))
    log_cpm <- do.call("cbind", log_cpm)
    colnames(log_cpm) <- unique(group)
                  
    # Compute log2fc between contrast and maximum group 
    max_log2fc <- data.frame(max_log2fc = log2(cpm[, contrast] / apply(cpm[, colnames(cpm) != contrast], 1, max)))
    rownames(max_log2fc) <- rownames(cpm)
    
    # filter log_cpm and max_log2fc by deg top hits
    log_cpm <- log_cpm[rownames(log_cpm) %in% rownames(deg), ]
    max_log2fc <- max_log2fc[rownames(max_log2fc) %in% rownames(deg), , drop = FALSE]
                      
    # Sort log_cpm by max_log2fc 
    log_cpm <- log_cpm[arrange(max_log2fc, desc(max_log2fc)) %>% rownames(), ]
                      
    # Filter for min expression in contrast group 
    log_cpm <- log_cpm[log_cpm[, contrast] > min_contrast_log2cpm, ]

    # Filter log_cpm top hits 
    if (nrow(log_cpm) > top) {log_cpm <- log_cpm[1:top, ]}
                      
    # ComplexHeatmap color ramp
    color_range <- max(abs(log_cpm))
    color_break <- seq(0, color_range, 0.01)
    color_ramp <- viridis(length(color_break), option = "magma")

                  
    # Heat map 
    hm <- grid.grabExpr(draw(ComplexHeatmap::pheatmap(
        mat           = as.matrix(log_cpm),
        main          = contrast,
        fontsize_row  =  10,
        scale         = "none",
        cluster_rows  = TRUE,
        cluster_cols  = FALSE,
        cellwidth      = 12, 
        cellheight     = 12, 
        clustering_distance_rows = "euclidean",
        clustering_distance_cols = "euclidean",
        clustering_method        = "complete",
        show_row_dend = FALSE, 
        show_rownames = TRUE,
        show_colnames = TRUE,
        color         = color_ramp,
        breaks        = color_break, 
        border_color  = NA)))

    
    return(hm)
    
}

# Contrast matrix for one vs average of all other explanatory variables 
contrast_one_vs_average <- function(design) {
    
    exp_variables <- colnames(design) 
    contrast <- list()
    for(i in seq_along(exp_variables)) {
        design_string <- paste(exp_variables[i], "-", "(", paste(exp_variables[-i], collapse = "+"), ")/", length(exp_variables[-i]), sep = "")
        contrast[[i]] <- makeContrasts(contrasts = design_string, levels = exp_variables)
    }
    
    return(contrast)
    
}

# Compute DEG matrix needs a contrast, design and count matrix object
deg <- function(
    
    cnt, design, contrast, 
    filter_min_cpm = 0.5,
    filter_min_expressed_samples = 2, 
    normalization_method = "TMM",
    adj_p_cutoff = 0.05
    
) {
    
    # Filter and sort cnt by design matrix
    cnt <- cnt[, colnames(cnt) %in% rownames(design)]
    cnt <- cnt[, match(rownames(design), colnames(cnt))]
    
    # Create a DGEList object from cnt and normalize
    cnt <- DGEList(cnt)
    cnt <- calcNormFactors(cnt, method = normalization_method)
    
    # Keep only rows where filter_min_expressed_samples or more samples out of group1 and group2 must be greater than filter_min_cpm
    cnt <- cnt[rowSums(edgeR::cpm(cnt, normalized.lib.sizes = TRUE, log = FALSE, prior.count = 0) > filter_min_cpm) >= filter_min_expressed_samples, ] # cpm default is normalized.lib.sizes = TRUE and log = FALSE
    
    # Limma fit to fit explanatory variable
    fit <- lmFit(voom(cnt, design = design, normalize.method = "none"), design)
    # For a linear model fit, compute moderated t-statistics, moderated F-statistic, and log-odds of differential expression by empirical Bayes moderation of the standard errors towards a common value.
    eb <- eBayes(contrasts.fit(fit, contrast))
    
    # Fetch topTable
    m = topTable(eb, adjust = "fdr", number = Inf)   # This should fetch every gene/probe sorted by adj p
    m = m[m$adj.P.Val < adj_p_cutoff, ]

    # select column subset and remove dots in adj.P.Val colname
    m = m[c("logFC", "adj.P.Val","AveExpr")] # missing logFC
    colnames(m)[2] = "adjPValue"

    return(m)
}

Import data

cnt <- read.delim(cnt_file, row.names = 1)
meta <- read.delim(meta_file)

Select samples by meta

# Select sample by cell_name 
cell_name <- c(
    
    "Fob", "MZB", "GCB", "SplPlsC", #B-cells SPLEEN 
    "BasoBM",  #Basophil BM 
    "cDC2", "pDC", #DC conventional and plasmacytoid SPLEEN
    "EoP", "Eo", #Eosinophil progenitor and adult BM
    "MEP", "PreCFUE", "CFUE", "EryBlPB", "EryBlPO", "Retic", #Erythrocyte BM
    "MonoBM", "Mac", #Monocyte BM and macrophage PERITONEAL  
    "Mast", #Mast cell PERITONEAL 
    "MegTPO", #Megakaryocyte cultured from BM
    "STHSC", "LSK", "MPP", #Multi potential progenitor BM
    "NeutBM", #Neutrophil BM
    "NK", #Natural killer cell SPLEEN 
    "CMP", #Common myeloid Progenitor BM
    "GMP", #Granulocyte macrophage progenitor BM
    "CLP", #Common lymphoid progenitor BM
    "CD4T", "CD8T" #T-cells total SPLEEN and LYMPHNODES
    
)

meta <- meta[meta$cell_name %in% cell_name, ]

Select sample and order count matrix by meta

# Filter and order count matrix columns by meta data samples 
cnt <- cnt[, colnames(cnt) %in% meta$sample_id]
cnt <- cnt[, match(meta$sample_id, colnames(cnt))]

Reset meta sample id and cnt column names

meta$sample_id <- make.names(paste0(meta$fine_labels, ".", seq(1:nrow(meta))))
colnames(cnt) <- meta$sample_id

Translate count matrix ENSEMBL to Symbol

# Get ensembl to symbol relationship from biomaRt 
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
ens_to_sym <- getBM(filters = "ensembl_gene_id", attributes = c("ensembl_gene_id", "mgi_symbol"), values = rownames(cnt), mart = mart)

# Remove duplicates 
ens_to_sym <- ens_to_sym[!ens_to_sym$ensembl_gene_id %>% duplicated, ]
ens_to_sym <- ens_to_sym[!ens_to_sym$mgi_symbol %>% duplicated, ]

# Filter data by ensembl id in reference 
cnt <- cnt[rownames(cnt) %in% ens_to_sym$ensembl_gene_id, ]

# order ens_to_sym ensembl id by rownames of data 
ens_to_sym <- ens_to_sym[match(ens_to_sym$ensembl_gene_id, rownames(cnt)), ]

# Set mgi_symbol as data rownames
rownames(cnt) <- ens_to_sym$mgi_symbol

Filter count matrix by genes

# Keep only genes with counts greater 0
cnt <- cnt[rowSums(cnt) > 0, ]

Set factor level

# Set factor level for labels 
meta$main_labels <- factor(meta$main_labels, levels = names(color$main_labels_haemosphere))
meta$fine_labels <- factor(meta$fine_labels, levels = names(color$fine_labels_haemosphere))

PCA on all data

options(repr.plot.width = 10, repr.plot.height = 5)

# create a DGEList object from cnt and normalize
cnt_pca <- cnt[rowSums(edgeR::cpm(cnt, normalized.lib.sizes = FALSE, log = FALSE) > 0.5) >= 2, ]
pca <- prcomp(t(edgeR::cpm(cnt_pca, normalized.lib.sizes = FALSE, log = TRUE)), scale = TRUE, center = TRUE)
pca_1 <- pca_scattern(pca, meta$main_labels, color$main_labels_haemosphere)
pca_2 <- pca_scattern(pca, meta$fine_labels, color$fine_labels_haemosphere)

pca_1 + pca_2

Design and contrast matrix

Explanatory variable (factor): Cell lineage (design_cl) or cell type (design_ct)
Indipendent meassurements (character): Cell type with batch information
Statistical model: Means model to get all possible combinations of levels
Design matrix: Coded as model.matrix(~ 0+explanatory_variable) which is here equivalent to model.matrix(~ explanatory_variable)

Compare main label

meta_ml <- meta
meta_ml$main_labels <- make.names(meta_ml$main_labels)

# Make design for means model parameterization
design_ml = model.matrix(~0+as.factor(meta_ml$main_labels))
colnames(design_ml) = unique(meta_ml$main_labels)
rownames(design_ml) = meta_ml$sample_id

# Create contrast 
contrast_ml <- contrast_one_vs_average(design_ml)

# Compute deg for all contrast combination 
deg_ml <- list()
for(i in seq_along(contrast_ml)) {deg_ml[[i]] <- deg(cnt = cnt, design = design_ml, contrast = contrast_ml[[i]])}

options(repr.plot.width = 20, repr.plot.height = 40)

# Plot heat map of marker genes
hm_ml <- lapply(seq_along(unique(meta_ml$main_labels)), function(i) top_hm(deg = deg_ml[[i]], cnt = cnt, meta = meta_ml, level = "main_labels", contrast = unique(meta_ml$main_labels)[i], top = 50, min_contrast_log2cpm = 2))
gridExtra::arrangeGrob(grobs = hm_ml, ncol = 4) %>% grid::grid.draw()

Compare fine labels

meta_fl <- meta
meta_fl$fine_labels <- make.names(meta_fl$fine_labels)

# Make design for means model parameterization
design_fl <- model.matrix(~0+as.factor(meta_fl$cell_name))
colnames(design_fl) <- unique(meta_fl$cell_name)
rownames(design_fl) <- meta_fl$sample_id

# Create contrast 
contrast_fl <- contrast_one_vs_average(design_fl)

# Compute deg for all contrast combination 
deg_fl <- list()
for(i in seq_along(contrast_fl)) {deg_fl[[i]] <- deg(cnt = cnt, design = design_fl, contrast = contrast_fl[[i]])}

options(repr.plot.width = 20, repr.plot.height = 90)

# Plot heat map of marker genes 
hm_fl <- lapply(seq_along(unique(meta_fl$fine_labels)), function(i) top_hm(deg = deg_fl[[i]], cnt = cnt, meta = meta_fl, level = "fine_labels", contrast = unique(meta_fl$fine_labels)[i], top = 50, min_contrast_log2cpm = 1))
gridExtra::arrangeGrob(grobs = hm_fl, ncol = 3) %>% grid::grid.draw()

Compare erythoblast stages

meta_ery <- meta[meta$main_labels == "Ery", ]
meta_ery$fine_labels <- make.names(meta_ery$fine_labels)

# Make design for means model parameterization
design_ery <- model.matrix(~0+as.factor(meta_ery$cell_name))
colnames(design_ery) <- unique(meta_ery$cell_name)
rownames(design_ery) <- meta_ery$sample_id

# Create contrast 
contrast_ery <- contrast_one_vs_average(design_ery)

# Compute deg for all contrast combination 
deg_ery <- list()
for(i in seq_along(contrast_ery)) {deg_ery[[i]] <- deg(cnt = cnt, design = design_ery, contrast = contrast_ery[[i]])}

options(repr.plot.width = 20, repr.plot.height = 10)

# Plot heat map of marker genes 
hm_ery <- lapply(seq_along(unique(meta_ery$fine_labels)), function(i) top_hm(deg = deg_ery[[i]], cnt = cnt, meta = meta_ery, level = "fine_labels", contrast = unique(meta_ery$fine_labels)[i], top = 50, min_contrast_log2cpm = 1))
gridExtra::arrangeGrob(grobs = hm_ery, ncol = 6) %>% grid::grid.draw()

Compare multipotent progenitor cells

meta_mpp <- meta[meta$main_labels == "MPP", ]
meta_mpp$fine_labels <- make.names(meta_mpp$fine_labels)

# Make design for means model parameterization
design_mpp <- model.matrix(~0+as.factor(meta_mpp$cell_name))
colnames(design_mpp) <- unique(meta_mpp$cell_name)
rownames(design_mpp) <- meta_mpp$sample_id

# Create contrast 
contrast_mpp <- contrast_one_vs_average(design_mpp)

# Compute deg for all contrast combination 
deg_mpp <- list()
for(i in seq_along(contrast_mpp)) {deg_mpp[[i]] <- deg(cnt = cnt, design = design_mpp, contrast = contrast_mpp[[i]])}

options(repr.plot.width = 20, repr.plot.height = 10)

# Plot heat map of marker genes 
hm_mpp <- lapply(seq_along(unique(meta_mpp$fine_labels)), function(i) top_hm(deg = deg_mpp[[i]], cnt = cnt, meta = meta_mpp, level = "fine_labels", contrast = unique(meta_mpp$fine_labels)[i], top = 50, min_contrast_log2cpm = 1))
gridExtra::arrangeGrob(grobs = hm_mpp, ncol = 6) %>% grid::grid.draw()

Compare restricted progenitor cells

meta_rpp <- meta[meta$main_labels == "RPP", ]
meta_rpp$fine_labels <- make.names(meta_rpp$fine_labels)

# Make design for means model parameterization
design_rpp <- model.matrix(~0+as.factor(meta_rpp$cell_name))
colnames(design_rpp) <- unique(meta_rpp$cell_name)
rownames(design_rpp) <- meta_rpp$sample_id

# Create contrast 
contrast_rpp <- contrast_one_vs_average(design_rpp)

# Compute deg for all contrast combination 
deg_rpp <- list()
for(i in seq_along(contrast_rpp)) {deg_rpp[[i]] <- deg(cnt = cnt, design = design_rpp, contrast = contrast_rpp[[i]])}

options(repr.plot.width = 20, repr.plot.height = 10)

# Plot heat map of marker genes 
hm_rpp <- lapply(seq_along(unique(meta_rpp$fine_labels)), function(i) top_hm(deg = deg_rpp[[i]], cnt = cnt, meta = meta_rpp, level = "fine_labels", contrast = unique(meta_rpp$fine_labels)[i], top = 50, min_contrast_log2cpm = 1))
gridExtra::arrangeGrob(grobs = hm_rpp, ncol = 6) %>% grid::grid.draw()

Create SingleR object

# Create SingleR object from Haemosphere RNA-seq data
ref <- SummarizedExperiment(list(counts = cnt, logcounts = edgeR::cpm(cnt, log = TRUE, prior.count = 1)))
ref$label.main <- as.character(meta$main_labels)
ref$label.fine <- as.character(meta$fine_labels)

Annotate SPLENO single cell data with Haemosphere SingleR object

# Import Seurat test data
so <- readRDS(so_file)

# Seurat to SingleCellExperiment 
sce <- SingleCellExperiment(list(counts = so@assays$RNA@counts))

# Predict labels
label_main <- SingleR::SingleR(test = sce, ref = ref, labels = ref$label.main, assay.type.test = "counts", de.method = "classic")
label_fine <- SingleR::SingleR(test = sce, ref = ref, labels = ref$label.fine, assay.type.test = "counts", de.method = "classic")

# Add labels to Seurat object
label_main_meta <- as.data.frame(label_main) %>% dplyr::select(pruned.labels, tuning.scores.first) %>% dplyr::rename(main_labels = pruned.labels, main_delta_score = tuning.scores.first)
label_fine_meta <- as.data.frame(label_fine) %>% dplyr::select(pruned.labels, tuning.scores.first) %>% dplyr::rename(fine_labels = pruned.labels, fine_delta_score = tuning.scores.first)

so <- AddMetaData(so, cbind(label_main_meta, label_fine_meta))

# Set factor level for labels 
so$main_labels <- factor(so$main_labels, levels = names(color$main_labels_haemosphere))
so$fine_labels <- factor(so$fine_labels, levels = names(color$fine_labels_haemosphere))

Save SingleR Haemosphere annotation

# write.csv(cbind(label_main_meta, label_fine_meta), "data/object/components/annotation/singler_haemosphre_counts.csv")
so <- AddMetaData(so, read.csv("data/object/components/annotation/singler_haemosphre_counts.csv", row.names = 1))

Plot SingleR delta score with label count

The main labels are not predicted by SingleR! We just use the top annoation for grouping e.g. the top lable for a particular cell might differ.

options(repr.plot.width = 15, repr.plot.height = 8)

hm <- function(data, fill) {
    
    set <- data
    
    # Mean delta score for fine labels per tissue group 
    data <- dplyr::group_by(data, fine_labels, tissue) %>% dplyr::summarise(mean_delta_score = mean(fine_delta_score))
    # Add main_label information for fine_labesl 
    data <- dplyr::left_join(meta[, c("fine_labels", "main_labels")] %>% unique(), data, by = "fine_labels")

    # Reshape to get scores by tissue 
    data <- reshape2::dcast(data, fine_labels+main_labels~tissue, value.var = "mean_delta_score")

    # Set NA scores to 0 for plotting 
    data$Myeloid[is.na(data$Myeloid)] <- 0
    data$Progenitor[is.na(data$Progenitor)] <- 0

    data$Myeloid[is.na(data$fine_labels)] <- 0
    data$Progenitor[is.na(data$fine_labels)] <- 0
    
    # Set factor levels
    data$main_labels <- factor(data$main_labels, levels = names(color$main_labels_haemosphere))
    data$fine_labels <- factor(data$fine_labels, levels = names(color$fine_labels_haemosphere))

    # Get matrix for heatmap 
    m <- t(as.matrix(data[, c("Myeloid", "Progenitor")]))
    colnames(m) <- data$fine_labels

    # Column annotation (not used)
    col_anno <- data.frame(data$main_labels)
    colnames(col_anno) <- c('Cell lineage')
    col_colors <- list("Cell lineage" = color$main_labels_haemosphere)
    col_legend <- list("Cell lineage" = FALSE)

    col_anno <- HeatmapAnnotation(
        which = "col", 
        df = col_anno, 
        col = col_colors, 
        show_legend = c(FALSE)
    )

    # Color ramp 
    color_ramp <- viridis(200, option = "magma")

    # Count histogram column annotation
    cell_name_count_p <- dplyr::filter(set, tissue == "Progenitor") %>% group_by(fine_labels) %>% dplyr::summarise(n = n()/nrow(.))
    cell_name_count_m <- dplyr::filter(set, tissue == "Myeloid") %>% group_by(fine_labels) %>% dplyr::summarise(n = n()/nrow(.))

    cell_name_count_p <- dplyr::left_join(data, cell_name_count_p, by = "fine_labels")
    cell_name_count_m <- dplyr::left_join(data, cell_name_count_m, by = "fine_labels")

    col_bar = HeatmapAnnotation(
        "Myeloid fraction" = anno_barplot(cell_name_count_m$n, gp = gpar(border = NA, fill = fill, lty = "blank")),
        "Progenitor fraction" = anno_barplot(cell_name_count_p$n, gp=gpar(border = NA, fill = fill, lty = "blank"))
    )

    hm <- grid.grabExpr(draw(ComplexHeatmap::pheatmap(
            mat            = m,
            name           = "Mean delta score",
            main           = "",
            fontsize_row   =  12,
            scale          = "none",
            cellwidth      = 20, 
            cellheight     = 20, 
            cluster_rows   = FALSE,
            cluster_cols   = FALSE,
            show_rownames  = TRUE,
            show_colnames  = TRUE,
            color          = color_ramp,
            column_split   = data$main_labels,  
            top_annotation = col_bar, 
            column_title_rot = 90,
            border_color   = NA)))

    return(hm)
    
    
}

# Fine label delta score heatmap with fraction of cells histogram annotation 
hm_nacl <- hm(subset(so, subset = treatment == "NaCl")@meta.data, color$treatment["NaCl"])
hm_cpg <- hm(subset(so, subset = treatment == "CpG")@meta.data, color$treatment["CpG"])

gridExtra::arrangeGrob(grobs = list(hm_nacl, hm_cpg), nrow = 2, ncol = 1) %>% grid::grid.draw()
ggsave(hm_nacl, filename = "result/plot/progenitor_nacl_haemosphere_deltascore.png", width = 10, height = 5)
ggsave(hm_cpg, filename = "result/plot/progenitor_cpg_haemosphere_deltascore.png", width = 10, height = 5)

Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from rna_umap_ to rnaumap_”
`summarise()` has grouped output by 'fine_labels'. You can override using the `.groups` argument.

Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from rna_umap_ to rnaumap_”
`summarise()` has grouped output by 'fine_labels'. You can override using the `.groups` argument.

Plot UMAP with SingleR Haemosphere annotation

reduction <- "rna_umap_nno"

options(repr.plot.width = 10, repr.plot.height = 15)

dplot_1 <- DimPlot(so, reduction = reduction, group.by = "rna_snn_res.0.8", label = TRUE) & 
    ggtitle("Cluster") & xlab("UMAP 1") & ylab("UMAP 2") &
    theme(aspect.ratio = 1, legend.position = "none") & 
    guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))

dplot_2 <- DimPlot(so, reduction = reduction, group.by = "tissue", label = FALSE) & 
    ggtitle("Tissue") & xlab("UMAP 1") & ylab("UMAP 2") &
    theme(aspect.ratio = 1, legend.position = "bottom") & 
    scale_color_manual(values = color$tissue, na.value = "dark gray") & 
    guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))

dplot_3 <- DimPlot(subset(so, subset = tissue == "Myeloid"), reduction = reduction, group.by = "main_labels", label = FALSE) & 
    ggtitle("Myeloid") & xlab("UMAP 1") & ylab("UMAP 2") &
    theme(aspect.ratio = 1, legend.position = "bottom") & 
    scale_color_manual(values = color$main_labels_haemosphere, breaks = names(color$main_labels_haemosphere), na.value = "dark gray") & 
    guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))

dplot_4 <- DimPlot(subset(so, subset = tissue == "Progenitor"), reduction = reduction, group.by = "main_labels", label = FALSE) & 
    ggtitle("Progenitor") & xlab("UMAP 1") & ylab("UMAP 2") &
    theme(aspect.ratio = 1, legend.position = "bottom") & 
    scale_color_manual(values = color$main_labels_haemosphere, breaks = names(color$main_labels_haemosphere), na.value = "dark gray") & 
    guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))

dplot_5 <- DimPlot(subset(so, subset = tissue == "Myeloid"), reduction = reduction, group.by = "fine_labels", label = FALSE) & 
    ggtitle("Myeloid") & xlab("UMAP 1") & ylab("UMAP 2") &
    theme(aspect.ratio = 1, legend.position = "bottom") & 
    scale_color_manual(values = color$fine_labels_haemosphere, breaks = names(color$fine_labels_haemosphere), na.value = "dark gray") & 
    guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))

dplot_6 <- DimPlot(subset(so, subset = tissue == "Progenitor"), reduction = reduction, group.by = "fine_labels", label = FALSE) & 
    ggtitle("Progenitor") & xlab("UMAP 1") & ylab("UMAP 2") &
    theme(aspect.ratio = 1, legend.position = "bottom") & 
    scale_color_manual(values = color$fine_labels_haemosphere, breaks = names(color$fine_labels_haemosphere), na.value = "dark gray") & 
    guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))

dplot <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + plot_layout(ncol = 2)
dplot
ggsave(dplot, filename = "result/plot/dimplot_haemosphere.png", width = 10, height = 15)

Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from rna_umap_ to rnaumap_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from rna_umap_ to rnaumap_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from rna_umap_ to rnaumap_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from rna_umap_ to rnaumap_”